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In this paper a novel computational technique for finite discrete approximation of continuous dy- 
namical systems suitable for a significant class of biochemical dynamical systems is introduced. The 
method is parameterized in order to affect the imposed level of approximation provided that with 
increasing parameter value the approximation converges to the original continuous system. By em- 
t— 5 ploying this approximation technique, we present algorithms solving the reachability problem for 

^^ biochemical dynamical systems. The presented method and algorithms are evaluated on several ex- 

emplary biological models and on a real case study. This is a full version of the paper published in 
the proceedings of CompMod 201 1 . 

1 Introduction 

CO 

O 

Under the modern holistic paradigm provided by systems biology 0, genome-scale knowledge of in- 
dividual components is combined with knowledge of interactions underlying the physiology of living 
J> organisms. The central goal of systems biology is to integrate all available biological data and to re- 

construct executable models [19] which allow to investigate complicated behaviour emerging from the 
underlying biochemistry. An important dimension is the quantitative aspect of the data and processes 

l/~) being modeled. 

f — With respect to [18], we consider biological models to be captured by the notion of a biochemical 

dynamical system consisting of variables describing a certain quantity of the respective species in time 
(e.g., number of molecules or molar concentration). Variable values evolve in time with respect to rules 
modeling the effect of reactions. The space of all possible configurations of variable values is referred as 

• ^h the state space. 

There exist several modeling approaches that differ in abstraction employed for modeling of time, 
variable values, and molecular interaction effects. The most commonly used approach concerns systems 
of ordinary differential equations (ODE) [28] where both time and model variables are interpreted as 
continuous quantities. Effects of interactions are modeled in terms of continuous deterministic updates 
of variables. Variable values represent molar concentrations of the species. In general, the ODE ap- 
proach relies on many physical and chemical assumptions simplifying thermodynamic conditions under 
which particular biochemical phenomena can be modeled correctly [25]. It is important to note that even 
simple interactions such as second order reactions lead to non-linear ODEs. However, under certain as- 
sumptions, biological systems make specific subclasses of general non-linear dynamical systems. Such 
a specialization motivated development of specific analysis techniques |fT8l [71 01 [24J . 



*The work has been supported by the Grant Agency of Czech Republic grant No. 201/09/1389 and by the Czech ministry 
of education intent No. MSM0021622419. 

. © L. Brim et al. 

_ „ , ™„ , This work is licensed under the 

CompMod 2011 „ A . u . T . 

Creative Commons Attribution License. 



2 Reachability in Biochemical Dynamical Systems 

Nevertheless, dimensionality and complexity of biological models preclude satisfactory application 
of analysis methods implying that to explore the model dynamics the only practicable method is numeri- 
cal simulation. Since numerical simulation generates an approximate solution (a trajectory) starting from 
a single initial point in the continuous state space, the scope of such exploration is limited to the partic- 
ular trajectory only. This is sufficient for "local" analysis provided that initial conditions are precisely 
known. However, studied systems are typically under-determined in terms of uncertain quantitative pa- 
rameters and initial conditions. Therefore generalization of the exploration scope is necessary to reveal 
and understand the complicated emergent behaviour. An important example of a problem which cannot 
be effectively solved by local methods is global temporal property - the problem to decide whether a 
given dynamical phenomenon, e.g., oscillation or variables correlation, is globally present/absent for all 
considered initial conditions lTT6l l9l. 

In this paper we limit ourselves to a subclass of dynamical phenomena representing reachability of 
a given portion of the state space. Example of a global temporal property problem that belongs to this 
subclass is to identify minimal or maximal concentration of species reachable from a particular set of 
initial conditions. 

In general, the reachability problem is undecidable due to unboundedness and uncountability of the 
state space. However, since concentrations of species cannot expand infinitely, state spaces of biologi- 
cal systems dynamics can be considered bounded in most cases. Analysis can be therefore considered 
indirectly on suitable finite discrete approximations of continuous state spaces [24, 3]. 

For a significant class of biochemical dynamical systems determined by multi-affine vector fields 
(i.e., affine in each variable), there has been developed an over-approximative abstraction technique 
based on partitioning the continuous state space by a finite set of rectangles. Rectangles determine states 
of a rectangular transition system representing the finite discrete (over)approximation of the continuous 
state space Ifl4l . as shown in Figure [T^. The rectangular abstraction has been employed in 11241 for 
reachability analysis and further elaborated by model checking methods in |6l. The results show that the 
extent of spurious behaviour introduced by the abstraction is typically very high thus limiting satisfactory 
application of the method. The problem is based mainly on the fact that a transition between any two 
individual rectangles over- approximates the vector field on the border between the rectangles (a so- 
called facet, see Figure [I])) provided that the information regarding which trajectories starting in an 
entry facet evolve through a particular exit facet is abstracted out. This causes the rectangular transition 
system to generate many rectangle sequences in which there is no corresponding trajectory of the original 
continuous system embedded. Moreover, the extent of such spurious behaviour is not directly eliminated 
by increasing the partition density. 

When analysing approximate models as in systems biology, the need for precise results critically 
required in systems verification can be relaxed provided that a suitable approximation of the solution can 
be even more efficient to obtain useful results. Henceforth, in the field of complex systems exhaustive 
techniques are often combined with approximative methods thus making a certain shift in the way of 
applying formal methods ||29l [TTl ITOl . 

1.1 Our Contribution 

We present a new technique for discrete approximation of biochemical systems with dynamics given by a 
system of ODEs with multi-affine right hand side. Our discrete approximation is not an exact abstraction 
wrt the original continuous system, but rather an approximation that approaches exact reachability with 
decreasing approximation granularity. While still assuming the rectangular partition at the background, 
we employ a measure that enables local quantification of the amount of trajectories evolving on a rectan- 
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gle in a particular facet-to-facet direction. To this end, every rectangle is augmented with a local memory 
representing the information at which part (entry set) of the entry facet it has been entered. On each 
entry set, we identify focal subsets from which all trajectories lead to the same exit facet. In Figure [jj, 
there are two different states of a quantitative discrete approximation automaton (QDAA) depicted. Both 
states share the same rectangle [1, 1.5] x [1, 1.5] and they differ in entry sets (marked yellow). The upper 
state with entry set {1.5} x [1, 1.5] has only one focal set - all trajectories from its entry set exit the state 
through the facet [1, 1.5] x {1}. The second state with entry set [1, 1.5] x {1.5} has two focal subsets 
made by the green and the red part of the entry facet, respectively. 

Transitions from a state with given entry set have weights assigned to themselves. Consider a tran- 
sition from a state A to a state B. The transition exists if there is a part P of the entry set of A such that 
the trajectories of ODE solutions go from P to B. Weight of a transition from A to B corresponds to 
the (n — 1 -dimensional) volume of P divided by the volume of the entry set of A. In this manner, the 
measure reflects amounts of trajectories proceeding in a particular direction. Rectangle regions related 



by weighted transitions make the QDAA which is a discrete-time Markov chain. (See Theorem 3.2 and 
its proof.) 

From a computational viewpoint, the continuous volumes are finitely approximated by discretization 
on a uniform grid. Local numerical simulations are employed to identify the entry regions and focal 
subsets. The density of facet discretization grid is considered as the method parameter. Because of com- 
bining numerical simulation with rectangular abstraction, the resulting QDAA makes neither an over- 
nor an under- approximation of the original continuous system. Since for every sequence of states the ap- 
proximate volume measure converges to the continuous volume with increasing discretization parameter, 
the parameter indirectly affects the correspondence between the original continuous behaviour and its ap- 
proximation. This makes the method sufficient for approximating reachability in complex biochemical 
dynamical systems. 

In general, the following main contributions are brought by this paper. 

1. A novel computational technique for finite discrete approximation of multi-affine dynamical sys- 
tems by means of QDAA. 



2. Showing that QDAA converges to the original continuous system behaviour. (See Theorem 3.3 
and its proof.) 

3. A reachability algorithm for QDAA. 

4. Evaluation on elementary models and an E. Coli case study. 

Since the most common application of the considered systems class is the domain of biochemical 
dynamical systems modeled directly by rules of mass action kinetics ll23l . evaluation of the method and 
algorithms is realized on biological models fitting this framework. 

1.2 Related Work 

Discrete approximation methods are commonly used in continuous and hybrid systems analysis (see 
for an overview regarding reachability) to handle the uncountability of the state space. Direct methods 
work on the original system and rely on a successor operation iteratively computing the reachable set 
whereas indirect methods abstract from the continuous model by a finite structure for which the anal- 
ysis is simpler. Our method belongs to the latter class, since it uses numerical simulations and creates 
the abstraction automaton. Considering a fixed set of initial conditions, there is a certain overhead with 



Reachability in Biochemical Dynamical Systems 



t i i 4 4 

i i i i 
i t i i i 


i l HI l 
1 J J t J 
j i h 1 if t. 


1 41 O f/U 

/ r/fj /' 


•'•'/, 


/////V/ 


//? / 









/ 
/ 










\ / 






{, 


, ' 


b 


dX 
it 

(IV 


= -0.1 ■ X 
= -0.3. Y 




1/2 


\/ \ 



OUT 1 1/2 



Figure 1 : (a) Vector field of a linear system partitioned by thresholds, (b) the principle of rectangular 
abstraction, (c) and quantification of the extent of over-approximation in terms of transition weights. 
The dashed line inside the rectangle demonstrates the approximate border separating trajectories exiting 
through different facets. 



generating states of the automaton in comparison with simple numerical simulations. However, the ad- 
vantage of constructing the automaton is obtaining a global view of the dynamics. Moreover, in addition 
to rectangular abstraction, the automaton is augmented with weighted transitions which represent quan- 
titative information describing volumes of subsets of initial conditions belonging to attraction basins of 
different parts of the phase space. 

An indirect method based on rectangular abstraction automaton making the finite quotient of the 
continuous state space has been employed, e.g., in E4l [Tl l3l. In general, these methods rely on re- 
sults lTT4l l20l and are applicable to (piece-wise) affine or (piece-wise) multi-affine systems. Although 
not addressed formally in this paper, our technique can be considered as a refinement of [24]. However, 
we focus on obtaining satisfactory approximate results eliminating the extent of spurious behaviour com- 
ing from conservativeness of rectangular abstraction. Our technique can be employed for the recognition 
of spurious behaviour of the rectangular abstraction transition system. 

The technique presented in ll27ll employes timed automata for the finite quotient of a continuous 
system as an alternative to piece-wise linear approximations. Another indirect technique adapted to 
multi-affine biological models is Ifl5l . The approach also employes rectangular abstraction, but results 
in less conservative reachable sets by means of polyhedral operations. In |f2j|8l there are techniques pro- 
posed for rectangular refinement that go towards reduction of over-conservativeness. These techniques 
work fine for linear systems while leaving the non-linear systems as a challenge. 

Direct methods are mostly based on hybridization realized by partitioning the system state space into 
domains where the local continuous behaviour is linearized [12J. This method, in an improved form, has 
been applied to non-linear biochemical dynamical systems ifTTl . In general, direct methods give good 
results for low-dimensional systems and small initial sets. In comparison with indirect approaches, they 
are computationally harder. From this viewpoint, our approach lies between both extremes. 



2 Preliminaries 



2.1 Basic definitions and facts 

Let N denote the set of positive integers, No the set NU {0}, and M ( | the set of nonnegative real numbers. 
For n £ N, denote W the standard ^-dimensional Euclidean space with standard topology and Euclidean 
norm |-| : R" — > Mj}~. For an arbitrary function / we use the common notation dom(f) for the domain of 
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f. 

For every i G {1 , . . . ,n} assume a^b\ G R such that a,- < bj. Denote / = nf=i [#«>£»] an n-dimensional 
closed interval in W and vol(I) the n-dimensional volume of I defined as vol(I) = Yl1=i {hi ~ a d- Further 
denote Inter(I) the interior of I, defined as the cartesian product of open intervals fl/Li ifliibi). 

For any XCI" denote X*(X) the Lebesgue outer measure (on W) of the set X. Basically A*(X) 
is the minimal nonnegative real number such that whenever X can be covered by a sequence of closed 
intervals in R" the sum of volumes of these intervals is greater then or equal to X*(X). (For precise 
definitions see l30l .) Note that A„*(X) < oo for every bounded set X and X*(l) = vol (I) for every n- 
dimensional interval /. 

Let n > 2, i < n,c G R. We use R" _1 (c) to denote the hyper-plane R" _1 (c) = {(xi , . . . ,x n ) G R" | x { = 
c}. Denote jt, : R" — >• R w_1 the projection omitting the jth variable, 7T,-((xi, . . . ,x n )) = (x\, . . . , x r _i ,jc ( - + i , ...x n 
Let X C M" _1 (c). We extend the notion of the (n — 1) -dimensional Lebesgue outer measure to such sets 
X and denote X*_ X (X) the (n — 1) -dimensional Lebesgue outer measure of JT,(X). 

Let / : R" — > W 1 be a continuous function (an autonomous vector field). We say that 

x = f(x) (1) 

is an autonomous ODE system. An important property of autonomous systems is the fact that if y(t ) is a 
solution of (fTl) on an open interval (a,b), then y(t + to) is also a solution (defined on interval (a — to,b — 

to))- 

A function / : R" — > W satisfies the Lipschitz condition locally on R", if for every x G R" there exists 
an open set U C R", x GU and a constant L 6 R such that for every two points x\ ,X2 € £/ the inequality 
|/(*i)-/(*2)| <£• |*i ~*2l holds - 

Theorem 2.1 (Trajectories of solutions of an autonomous system) Let (|7| be an autonomous ODE sys- 
tem, where f is defined on R" and let f satisfy the Lipschitz condition locally on R". Let x be an inex- 
tendible solution of system HI. Then dom(x) is an open interval, and for every point a G R" there exists 
exactly one trajectory of an inextendible solution x(t) of system ([7]) coming through a. 

Theorem 2.2 (Continuous dependency on initial conditions) Let f : R" — > R" be continuous on an 
open set £C1" with the property that for every yo G E, the initial value problem x = f(x) , x(0) = yo has 
a unique solution y{t) = T7(f,yo) f 7 ] w a function of variables t,yo ). Let w±,wj G R such that (w±,wj) 
is the maximal interval of existence ofy{t) = T](f,yo)- 

Then the bounds w±,Wj are (lower, resp. upper semicontinuous) functions ofyo in E and T}(t,yo) is 
continuous on the set {(t,yo) \ yo G E,w±(yo) <t< wj(yo)} C R" +1 . 

We restrict ourselves to multi-affine autonomous systems. That is, systems of the form ([TJ, such 
that the vector field / is a multi-affine function, defined as a polynomial of variables x\ , . . . ,x n G R" of 



degree at most one in every variable. The assumptions of Theorems 2.1 and 2.2 (from B21I ) are satisfied 
for systems of this class, therefore the properties stated in the above theorems can be used for reasoning 
about autonomous systems with multi-affine vector fields. 

2.2 Biochemical dynamical system 



According to [18], by a biochemical dynamical system we understand a collection of n biochemi- 
cal species interacting in biochemical reactions. Species concentrations are represented by variables 
x\,...,x n attaining values from Rt. If the stoichiometric coefficients in reactions do not exceed one 



Reachability in Biochemical Dynamical Systems 



and the reaction dynamics respects the law of mass action kinetics 1231 , the dynamical system can be 
described by a multi-affine autonomous system in the form ([TJ). 

In a biochemical dynamical system we are typically interested in a bounded part (^-dimensional in- 
terval) of the phase space in W. Further, we consider the phase space partitioned by a (non-uniform) 
rectangular grid. In particular, for each variable there is defined a finite set of thresholds, making the sys- 
tem partition. Thresholds determine (n — 1) -dimensional hyper-planes in R" and can be freely specified 
according to particular questions that should be addressed by the model analysis, e.g., specification of 
unsafe or attracting sets. Cells laid out by In adjacent threshold hyper-planes (cells are again intervals in 
W) are called hyper-rectangles, for short we refer to them as rectangles. 

Definition 2.1 Define a biochemical dynamical system (biochemical system for short) as a tuple S$ = 
(n,f,3?,J?c), where 

• n € N is the dimension of ' S$, 

• / : M" — > W is the multi-affine vector field of 'S3, 

ST = (7i, . . . ,T n ) is the partition of 33 where each Tj is a finite subset ojTKq, and define the set of 
rectangles given by 2? as 



• 



Rect(ST) = {Y[lj | Vj3a,b £ Tj : Ij = [a,b],Vc £Tj : c <a\J c >b}, 

7=1 

• <?c Q Rect(3 r ) is the set of initial conditions (initial set) of S3. 

Definition 2.2 Let S3 = {n,f, 2? ,^c) be a biochemical system and let H € Rect{3 r ) be a rectangle such 
that H = I\X ...Xl n , where /, = [ai,bj\. For every i E {1, . . . ,«} define the lower (resp. upper) facet of 
H wrt the ith variable: 

Facetj-(H) = {(xi,. . . ,x n ) GH\x{ = a,}, 

Facet] (H) = {(x\, . . . ,x n ) G H \ xt = bf\. 

Denote Facet Sj(H) the set of z'fh dimension facets of H, Facet st{H) = Facet f- (H) U Facet] (H), and 
Facets(H) the set of (all) facets of H, Facet s(H) = \J" =l FacetSi(H). 

Definition 2.3 LetH, H' £Rect(£?). We say that H is a neighbour of H', denoted H\x\H', if there exists 
F G Facets (H) such thatHDH' = F. 

3 Quantitative Discrete Approximation 

Given a biochemical system SB = (n,f, 3?,^P), we aim to define a finite automaton reflecting the be- 
haviour of S$, and for each state, to assign every transition a weight quantifying probability of proceeding 
to a particular successor. 

A state is defined as a pair (H,E) - a rectangle H, and a subset £ of a particular facet of H. The set 
E represents a so-called entry set, a region through which trajectories of the system ([TJ enter the interior 
of H. Intuitively, we can say that E encodes the history of previous evolution of the system from initial 
set J'c to H. Entry sets are either subsets of (n — 1) -dimensional facets of H or (in case of initial states) 
the whole ^-dimensional rectangle H. 

Since entry sets can be arbitrary sets in Euclidean space, we approximate them by a finite discrete 
structure. Each facet is provided with a uniform grid on which we approximate any subset of the facet 
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Figure 2: Example of a biochemical system with two species and two reactions. Dynamics given by 
a system of two ODEs and the system of thresholds are in the left part of the figure. Vector field is 
visualized in the middle, and its Rectangular Abstraction Transition System on the right. 

by the set of rectangular fragments, so-called tiles (Figure [5]). The grid is ra-dimensional or (n— 1) 
dimensional depending on the dimension of approximated entry sets. When following the trajectories of 
solutions of differential equations of the models dynamics in time, entry sets are identified by trajectories 
of solutions passing through them on their way from preceding rectangles. In following definitions we 
treat this intuitive perception of entry sets formally. 

Let K G N, let & = (n,f, ST, ^c) be a biochemical system, H G Rect(3"), and F G Facets{H) for all 
definitions and theorems from this section. 

Definition 3.1 LetHbeoftheformH = Y\" =l lj,where\/j:lj = [aj } bj\. Let B G {H}UFacets(H). Set 
either n 1 = n, ifB = H, or n! = (n— 1), ifB G FacetSj(H) for some 1 <i <n (in this case 3c G {a,-, ft,} : 

5cm? _1 (c);. 

Define the set of K-tiles of B as Tiles* (B) = {A C B \ A = Y\ n j=[ Aj}, where A; = {c}, if B G 
FacetSi[H), and otherwise ( j ^ i or B = H ) A y is a closed interval in M^ oftheform [a y + -^{b y — a y) ,a y + 
-^— - (bj —ay)], where for all j G {l,...,n},j ^i the nonnegative integer kj G No satisfies kj < K. 

The following definition introduces the notion of general entry sets. 

Definition 3.2 Define the set of entry points into a rectangle H through facet F, as the set Entry '(F,H) = 
{ y Q G F | 3 a trajectory y(t) of a solution of Q such that y(0) = y and 3e > : y(t) G H for \/t G (0, e) } . 

Next we define the approximation of entry sets on a grid of K- tiles. Additionally, we define the 
respective (discrete) volume measure of a set (see Figure |3]c),d)). 

Definition 3.3 Let X c H. Let n' = n — 1, if there exists i G {1, . . -,n},F G Facet S{(H) such that X C F, 
and let n' = n, otherwise. Let M = F, ifX C F, and let M = H, if there is no such facet F. Define 

• the set of K"-tiles approximating the set X as 

Tiles*, (X) = I A G Tiles*, (M) \ K ^* ] > ~ j, 



• the rectangular fc-grid measure of the setX as X*,{X) = Y,AeTUes K ,(x) vol(A). 

The following definition declares the set of all discretized entry sets for a given rectangle. 
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Figure 3: a) Let H = [0,2.5] x [0,2.5] be a rectangle. The blue areas depict elements of Tiles\(H). 

b) Let F = FacetJ(H) = {2.5} x [0,2.5]. The red line segments are elements of Tiles\{F). 

The set EntrySets?, (H ) has 2 +4 • (1 + Q) + (\) ) = 30 elements: 0, H, and 7 for every facet of H (the facet 
itself, 3 segments and 3 unions of pairs of segments of the facet). 

c) Let X be a subset of H (the shaded polygon). Let K = 5. 

d) The set of K— tiles approximating X is the set of five blue intervals (each satisfying that at least half of 
its area is in X). The cardinality of Tiles$(X) is 5. Thus X%(X) = 5 • (0.5 • 0.5) = 1.25. 



Definition 3.4 For H, define set of (approximate) entry sets EntrySets K (H) = 

iEQH\E = Q)\JE = H\j3F£ Facets(H),£Q Tiles K (Entry(F,H)) :E = \Jg\. 

For an example of a set of (approximate) entry sets of a rectangle see Figure [3] a),b). Note that set 
of approximate entry sets is always finite. Further note that also the empty set and the entire rectangle 
are considered as entry sets. These represent singular cases needed in the subsequent construction of 
the automaton. In particular, states with the empty entry set approximate fixed point behaviour not 
leaving the rectangle (steady state memory) whereas the rectangle-form entry set is employed for initial 
rectangles. 

Definition 3.5 Let E G Entry Sets K (H),H' G Rect(£?),F' G Facets(H) such that H' m H,F' =HDH'. 

Define the focal subset of E on H targeting F', denoted Focal(H,E,F'), as the set of all yo G E 
such that there exist £,s',c > and a trajectory of a solution y(t) of system ([7]) with inital conditions 
y(0) = yo satisfying y{t) GHfort G (0,c),y(t) G Inter{H) for t G (0,e),y(c) £F', andy(t) G Inter (H') 
for t G (c, c + e'). Let ExitSet(H,E,F') denote the set of all such (targeted) points y(c) G F'. 

Define focal subset of E on H not leaving H, Focal(H,E,d>), as the set of all points yo £ E such 
that there exists a trajectory of a solution y(t) of system (|7|) with initial conditions y(0) = yo satisfying 
y{t) ^H for all t > 0. 

Next we define the successor function for any pair (H,E) and subsequently the quantitative discrete 
approximation automaton. 

Definition 3.6 Let E G EntrySets K (H). Define the successors of {H,E) as the set of pairs (H' ,E'} with 
H' G Rect(,^),E' G EntrySets K (H') such that 

Succs((H,E}) = {(H' , E 1 ) I H',E r satisfy one of conditions 1. — 3. below} 

1. H 1 [xi H, E / 0. Denote F 1 the facet of H satisfying F 1 = H n H 1 . Let n' = n, if E = H, and 
n' = (n — 1), otherwise. Moreover, E' = \JTiles K (ExitSet(H ,E ,F')) and \*, [Focal(H ,E ,0)) > 0. 

2. H' = H, E / 0, and E' = 0. Further, it holds that either E CF and X*_ x (Focal(H,E,Q))) > 0, or 
E = H and X,f (Focal (H,E,®)) > 

3. H' =H andE' =E = ®. 
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Definition 3.7 (The Quantitative Discrete Approximation Automaton) Let k, S3 be as above. The 
quantitative abstraction automaton QDAA K {38) of a biochemical system 38 with parameter K is a tu- 
ple QDAA K (SS) = (S,J?c,8,p), where 

• the set of states S = {{H,E) \ H <= Rect(^),E G EntrySets K (H)} , 

• the set of initial conditions Iq = {{H,H} \ H G J^c}* 

• the transition function 8 : S — > 2 s is defined as 8((H,E)) = Succs((H,E}), 

• the weight function p : S x S — > [0, 1] is defined by the following expression, where S = (H,E),S' = 
(H' ,E'). Suppose n' = n, in case E = H, and n 1 = n — \, otherwise. 



p(S,S') 



X*,(Focal(H,E,(d)) 



LAEFacets(H)U{Q>} K' { F OCal(H ,E , A)) 

k^(Focal{H,E,F')) 

HAeFacets(H)u{@} K' (Focal(H,E,A)) 

0, otherwise 



if H=H',E = E' =%, 
if H = H' } E^<d,E' = 0, 

if Htxi H'.E' C F' =HDH\ 



Example 3.1 Assume the biochemical system from Figure H\ See Figured a) for an example of fo- 
cal subsets described below. Let R = [0,2.5] x [2.5,5] be a rectangle and let Fq = FacetJ(R),Fi = 
FacetJ (R),F2 = Facets (R) , F3 = Facetf(R). For the state (R,Fq) the focal set ofF\ equals Fq, whereas 
Focal (Fq) = Focal (F2) = Focal (F^) = 0. 

Let H = [0,2.5] x [0,2.5] and F = FacetJ (R). For the state (H,H) the set Focal(F) is the blue area 
inside H and Focal(Q)) is the yellow area. All the solutions of the biochemical systems dynamics with 
initial conditions in Focal((Z>) approach the yellow line of fixed points and stay in H forever. All the 
solutions starting in the blue area leave H infinite time through F. 

In the right part of Figure^is the set of reachable states of the quantitative discrete approximation 
automaton (QDAA) obtained from the biochemical system described in Figure^with initial conditions 
J?c = {[0,2.5] x [0,2.5]}. 

Let H,R be the same as above. Let S = [2.5,5] x [0,2.5] and let J?c = {H}. The QDAA successor 
states of(H,H) are (H,Q>) (a selfloop state) and (S,E) (where E denotes the K-tiles approximation of the 
red segment in Facetf(S)). For K — > 00 the weights of these two transitions approach the area ratios of 
yellow and blue regions ofH respectively. The only successor of (H ,0) is (by definition) itself. The state 
(S,E) has one successor (5,0), since all the trajectories beginning in E approach the line of fixed points 
and stay inside S forever. 

Therefore the set of concentrations reachable from initial rectangle H is [0, 5] x [0, 2.5]. See the rect- 
angular abstraction transition system from FigureVjwhere the set reachable from H is [0,5] x [0,2.5] U 
[2.5,5] x [2.5,5], although there exists no trajectory of a solution of the biochemical systems dynamics 
that starts in H and reaches a point inside [2.5,5] x [2.5, 5] . 

On the other hand, if K is too small, some behaviours of the system are not reflected in QDAA, 
because the set of K-tiles corresponding to the entry set may be empty. With finer partition into K-tiles 
smaller entry sets can be captured and approximation of the biochemical system by a QDAA is more 
realistic. 



In the Theorem 3.1 we ensure correctness of using the Lebesque measure in Definition 3.7 We 
ensure that there is no non-zero volume entry set such that all trajectories from this set lead to a facet 
without entering the interior of a neighbouring rectangle. 



Following notations and lemmas provide the background that is used in the proof of Theorem 3.1 
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Figure 4: a) Focal sets examples, b) QDAA example. 



Notation 3.1 Let us denote B e (x) the open sphere with centre x G M" and radius e > (B e (x) = {y£ 
M n ||jt-y| <£}). 

The following theorems from mathematical analysis and differential topology (see for example EL 
ED) will be used in the proofs. 

Lemma 3.1 (Peano) Let f : M M — > W be continuous on open set E C|" such that f possesses contin- 
uous first order partial derivatives. Then for every yo E E, the initial value problem x = f(x),x(0) = vo 
has a unique solution y(t) = T}(t,yo) and the unique solution has continuous first order derivative wrt 
the variable t (is of class C l ) on its open domain of definition. 

Lemma 3.2 The following statements about zero measure sets hold. (A set has Lebesgue outer measure 
zero iff it has Lebesgue measure zero.) 

1. IfS C R" , X* (S) = and f : W — > W is a smooth map (with continuous first order derivatives), then 

K(f(s))=o. 

2. Being of zero measure is a diffeomorphism-invariant property of subsets ofW. 

Lemma 3.3 (A version of Fubini theorem) Let U be a compact subset ofW 1 . Denote by Mf the subset 
{t} xr-'c W, and denote U t = UC\ Rf. If every set U, satisfies K-\( u t) = 0. then K ( u ) = °- 



Lemma 3.4 Let g 

hold. 



I" be a multiaffine function. Let F £ FacetSj(H). The following statements 



1. If there exists a set U C H such that X*(U) > and \/x £ U : g(x) = 0, then g = onH (and on W 1 ). 

2. If there exists a setV CF such that A*_j (V) > and Vx € V : g(x) = 0, then g = on F (and on the 
hypherplane M" (c) containing F). 

Proof The proof of statement 1. will be done by mathematical induction wrt n. 

For n = 1 is H a line segment, and the function g is a linear function of one variable. Either there is 
one point (a set of measure zero in M 1 ) where the function attains zero value, or the function is zero on 
the whole line. 
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For n = 2, H is a rectangle in plane. The multi-affine function of two variables is either zero on whole 
plane, or on two intersecting lines, or a hyperbolic curve or empty set. Therefore either the function is 
zero on the whole plane, or on a subset of plane of measure zero. Only the second case is compatible 
with g being nonzero on a set of non-zero Lebesgue outer measure. Assume that the statement holds for 
1 , 2, ... ,n — 1 and let us prove it for n. By Theorem [33] 

Proof of 2. is analogous, we can identify R"(c) with M" _1 and assume that g is a multiaffine function 
of n — 1 variables (with xt = c constant). □ 

Definition 3.8 Let U C H, let I be an (bounded or unbounded) interval in R. Define the /-trajectories 
unoin of U as the set 

Ei(U){x(s)\s €l,x(t) is a solution of system (fTl) defined on l,x(0) S £/}. 

Lemma 3.5 Let U C H be a (n — I) -dimensional closed disc in H (i.e. 

3i £ {1, . . .,n}3c G R38 > 03x £ H :U = {y £ R?(c)| \x- y\ < 8}, 

we do not assume there is a facet F ofH with U C F). 

Let X*_ l (U)>0 and let either My £U : f(y) ■e i >0 or\/y £U : f(y) ■ e t < 0. 

Then either 3U' C U with X*_y(U') > such that all the solutions of system nil stay forever in H, or 
there exists F' G Facets(H) such that X*_ 1 (F' n Sj . +oo ) (£/)) > 0- 

Moreover, in the second case, X*_ l (F' DW) > 0, where W is the set of points of the boundary ofH 
where the trajectories leave the rectangle H for the first time (a subset of the connected component of 
E[o )+00 )(E/) C\H containing U). 

Proof Consider the case, when there is no U' C U with X*_ 1 (U') > such that all the solutions of 
system ([TJ stay forever in H. 

Since all compactly supported vector fields are complete, we can restrict ourselves on the multi-affine 
vector field / defined on a compact neighbourhood of H and assume that flow of this (complete) vector 
field defines a one-parametric group of diffeomorphisms Fl{ : (t,yo) — > y(t), where y(t) is the respective 
solution of system [T] 

We consider a diffeomorphism cpu between Sr 0+oo )([/) C\H and a special subset V of U x M.. We 
prove that for U satisfying the assumptions of Lemma [33] the set V is of nonzero measure. Therefore its 



image by diffeomorphism q> must be of nonzero measure (due to Theorem 3.2 1. If on the other hand the 
set Sf-MjO] {F 1 H W) n H was of measure zero, its image in (pp'nw must be of zero measure also, but since 
E[o,+~) (U) C\H is the union of finitely many such sets Z/_„ i (F' C\W) HH, there has to exist F' such that 
the set E/^o] (F' n W) C\H has nonzero measure. □ 

Lemma 3.6 Let E G EntrySets K (F,H),E / 0. 

Then either 3U' C E with X*_i(U') > such that all the solutions of system nil stay forever in H, or 
there exists F' € Facets(H) such that 

X;_ y (F' r\^ +oo) (E^Entry(F,H))) > 0. 

Moreover, in the second case, X*_ l (F' D W) > 0, where W is the set of points of the boundary ofH 
where the trajectories leave the rectangle H for the first time (a subset of the connected component of 

Em +00 \ (ED Entry (F,H))DH containing the set EDEntry(F,H)). 
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Proof Consider again the case, when there is no U' C U with X*_ t (£/') > such that all the solutions of 
system ([TJ stay forever in H. 

The situation is easier if there exists a point x G E with f(x) ■ Vh{F) < 0. The multi-affine function / 
is continuous, therefore there exists am open sphere containing x such that f(y) ■ V#(F) < for all points 
y from this sphere. As a subset of the intersection of this sphere and F there exists a closed disc like U 
from the assumption of Lemma [33] and the statement of this lemma holds. 

Otherwise the whole set EnEntry(F,H) is a subset of F with nonzero (n — 1) -dimensional outer 



Lebesgue measure on which / = 0. By Lemma 3.4 / = on the whole (n — 1) -dimensional hypher- 
plane containing F. 

The idea of proof for this configuration is to use the defmitorial property of Entry(F,H), that insures 
existence of e z and a point y G Inter(H) reachable by a trajectory from a point z in EnEntry(F,H) in 
time e z with f(z) • V# (F) < (that is implied by the properties of multi-affine function /). By continuous 
dependency on the initial conditions we find a disc around z that cosists of points reachable from E n 



Entry (F,H) (continuous dependence) and satisfies assumptions of Lemma 3.5 Then we use Lemma 3.5 
and this lemma is proved. □ 

Theorem 3.1 Let E G Entry Sets K (H),E ^ 0. Further, let n' = n, ifE = H, and n' = (n — 1), otherwise. 
Then 

£ V n ,(Focal(H,E,A))>0, (2) 

AeFacets(H)U{®} 

Proof The proof will be divided into two parts. 
First, for E C F,E / 0. Second, for E = H. 



Proof of 1. The statement of theorem is obtained by using Lemma 3.6 

Proof of 2. Let E = H. In the case when all the trajectories with initial points in H stay forever in H 
the inequality X* [Focal(H ,E ,0)) > holds. 

In the case when there exists a point xq G H and a trajectory of a solution x(t) of system (fTl) with 
jc(0) =x ,3c,£ >0: V? € [0,c],x(t) €H,Vt £ (c,c + e)x(t) £ H. Let us denote y the point x(c+ f) £H. 



From Theorem 2.2 (continuous dependency on initial conditions) there exists an «-dimensional sphere 
B with centre xq such that trajectories of all solutions of ([[]) with initial point in the sphere B leave 
the rectangle H and continue into a neighbourhood of point y. The intersection HHB surely contains 



a disc satisfying assumptions of Lemma 3.5 and therefore there exists F' € Facets(H) such that the 
intersection of S[o i+ oo)(b) and F' is of nonzero measure, i.e. it has a subset of nonzero (n— 1) -dimensional 
measure in the interior of the facet F'. Then there exists an open sphere B' contained in B such that all 
trajectories of solutions with initial points in B' leave the rectangle H through the interior of facet F', and 

X*{Focal{H,E,F')) > X*(B') > 0. □ 

Theorem 3.2 The quantitative abstraction automaton QDAA K (^) of a biochemical system 33 is a dis- 
crete time Markov chain. 

Proof The number of states is finite, bounded by \Rect(&) \ ( 2 + 2n ■ (2 K " —V 

Sum of probabilities of transitions from one state (H,E) equals 1 for E = and AeFa "" ( "') — — e 

T.AeFacets(H)u{Qi} \* \F ocal(H ,E fi) J 



The later sum equals 1 , whenever its denominator is nonzero (it is the case because of Theorem 3.1). The 



probabilities of transitions from a given state are independent of previous states of the automaton. □ 

Finally, we provide a theorem suggesting that for sufficiently large values of parameter K, the rectan- 
gular K"-grid measure of a bounded set X contained in the phase space of biochemical system approaches 
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its Lebesque outer measure (Theorem 3.3 1. Following lemmas and definition will be used in the proof of 
this theorem. 

Let 38 = («,/, T,Ic) be a biochemical system, H G Rect(3?),F G Facets(H), and K G N throughout 
this section. 

Definition 3.9 Let X C ]T/ = i [ m i n (^') > max (7})] • Define X£(X) as the sum 






Analogously define X*_ l (X) as the sum 



£ ^(xnF). 

f eU«« c «(.y) K!C «"'( w ) 
Lemma 3.7 Lef m £{n — l,n}, let J be an m-dimensional interval in n-dimensional space, /C]T/ =1 [min(7}),max(7})]. 

limA^(7)=A; 7 (7) 

Proof For given K" there are at most 2rc • K"" -1 distinct tiles R satisfying R<^J and R n 7 7^ 0. 
The difference between A,£(/) and A*(/) = vo/(7) is bounded by 

l -£vol(R i )= l -2n-K»- [ vol(R 1 ) = 

= nK — = — 
K n K 

The expression approaches zero as K — > 00. □ 

Lemma 3.8 Let M be a positive integer number (M < 00). Let U = \J j=l l\ be a union ofM n-dimensional 
bounded rectangles, subsets o/TT)'=i [min(7}),max(7})]. Let this the intersection of every two intervals be 
of zero Lebesgue outer measure. Then 

Proof There are maximally M • 2nK n ~ l distinct tiles R satisfying R^-U and RC\ U 7^ 0. 
The difference between X£(U) and X^(U) = vol(U) is bounded by 

M~Y,vol{Ri) =M-2n-K"- 1 vol(Ri) = 

, V nV 

= MnK"- 1 —=M— 
K" K 

The expression approaches zero as K — > 00. □ 

Lemma 3.9 Let X be a (bounded) subset ofM." and let J\ ,J2, ■ ■ ■ be a (countable) sequence of intervals 
such that X C (JT=i ^j- Then there exists a sequence 1\ ,12, ■ ■ ■ of intervals with X* (/,• D/y) = Ofor each 
pair i 7^ j such that (J^=i h = U7=i //• 
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Proof The sequence of pairs of intervals from the first sequence is also countable. 

Replace every two rectangles that overlap by finitely many new non-overlapping rectangles. □ 

Lemma 3.10 Let X C fl/Li [ rmn (^)) max (^')] be a bounded subset ofW 1 with lebesgue outer measure 
\*(X) = reR+ 

Let e > 0. Let J\ : J%,...be a cover ofX by intervals such that YlfLi vol(Ji) — r<£. 

Let 8 > be a positive real number. Then there exists a number kg such that X^ij vol{Ji) £ (r — 
S,r + S). 

(That means TT=k 5 vol(Ji) < 8.) 

Proof The sequence of partiall sums of the sequence of volumes (vol(Ji))°° =l converges to r. That implies 
the the existence of such kg. □ 

Corollary 3.1 Let X and 8 be the same as in the previous lemma. Let Ii,h, ■ ■ ■ be a cover of X by 
countably many intervals with X*(IjDlj) = Ofor each pair i ^ j. Let kg be as in the previous lemma. 
Then there exist two subsets ofX denoted X' andX" such that X=X' UX", X* (X' f]X") = 0, X" C |jfii h 
and X' CUr=^ 4 

Proof The existence of X' and X" is obvious (can be defined as X" = (J*^ I t nX and X' = |J~ h hC\X 
respectively), there remains the proof of the equality X*(X' nX") = 0. 

The set X' C\X" contains only countably many intersections of intervals h,h, ■ ■ •> and these intersec- 
tion are of measure zero, therefore their union has also measure zero, X* (X' (IX") = 0. □ 

Lemma 3.11 Let X be a subset o/TTJLi [min(7]) , max(7})] with Lebesgue outer measure X* (X) = r. Then 
for every K the following statement is true: 

K(X) < 2r 

Proof For any tile J to be counted into the fc-grid measure of a set, the inequalities X* (X n J) > 
\vol(J) 44> 2X* (X n J) > vol(J) must hold. Therefore the set X with Lebesgue outer measure r can 
saturate at most a set of tiles of the overall volume 2r. □ 

Theorem 3.3 LetX CH<E Rect(^) (or more generally X C ]T/ = i [min(7}),max(7])], where (T u ...,T n ) = 
^ is the partition of 33). Then 

hmX n K (X)=X:(X). (3) 

Proof The proof uses the definition of outer Lebesgue measure and our aim is to prove that for every 
e > there exists such K that the difference \X*(X) - X*(X)\ <£. 

Let £ > be a positive real number, denote r = X* (X). 

For X in the type of interval the proof is easy (Lemma [3/7] ) . 

For general X consider sufficiently accurate (with less than | difference of sum of volumes and r) 
cover of X with countable collection of intervals whose interiors do not intersect (Lemma [379] >. 

There exist two subsets of X denoted X' and X" such that X = X' UX", X*(X' nX") = 0, such that 
the collection of intervals can be divided into a finite part I\, . . . ,4 and the remainder 4 + i, . . . such that 



the sum of volumes of the remainder is sufficiently small (less than |) (Corollary 3.1 1. 

For a bounded set Y with X* (Y) = s the inequality X*{Y) < 2s holds for every K" (Lemma [3. 11 [ >. 

For the finite set of k intervals we can find such K that | Y!j=i ^nifi) ~ Lf=i ^n (4)1 — I ( a PPlyi n g 
Lemma [377] finitely many times, taking the maximal of obtained values of K") and the overall difference 
of \ K (X) and X*(X) is less than £. □ 
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Note that the result applies also to the case with XCF£ Facets(H) and X*_ i , A*_ { . 

4 Algorithm 

This section introduces procedures for obtaining the reachable state space of the quantitative discrete 
approximation automaton. Algorithm [T] is a procedure of computing the set of reachable states. Al- 
gorithm [2] describes the computation of transitions from one state (i.e. successors) together with their 
weights using numerical simulations. 

The procedure of computing reachable state space (Algorithm [T]) is based on breadth first search. 
States corresponding to initial conditions of the biological system are enqueued first and a list of states 
already visited is maintained. The computation is always finite, because there are only finitely many 
possible states of the automaton and each of them can be at most once added and after the computation 
of its successors removed from the queue. 

Algorithm 1 Computing the set of reachable states 

Require: SB = («,/, &~, J'c), k e N 

Ensure: Reachable = set of all reachable states of the automaton QDAA K (SS) 



1: 


Reachable <- 


2: 


for all H e y c do 


3: 


s <- (H,H) 


4: 


Reachable <— Reachable U{s} 


5: 


Queue. pushBack(v) 


6: 


while Queue ^ do 


7: 


s <— Queue. firstElement 


8: 


A <— getSuccessors(.v) 


9: 


for all a e A do 


10: 


if a <£ Reachable then 


11: 


Reachable <— Reachable U{a} 


12: 


Queue. pushBack(a) 



13: return Reachable 



Computation of the successors (Algorithm |2|) of one state requires determining the rectangles and 
the entry sets of the successors and weights of the transitions. This can be done approximately using 
numerical simulations. We sample the entry set of the state and perform numerical simulations with 
the sampled points as initial conditions and the dynamics of the given biological system as the vector 
field. For each simulated trajectory we watch whether it leaves the rectangle before given maximal time 
interval elapses. If this is the case then the location of the exit point through which the trajectory leaves 
the rectangle is of interest. 

Entry sets of the successor states are also determined within Algorithm [2] If the successor is a 
selfloop state the entry set is empty. For a neighbouring rectangle successor with one common facet the 
entry set is computed using the exit points locations and more numerical simulations. From the set of 
exit points in a facet we can estimate the set of fc-tiles of the facet that surely have nonempty intersection 
with the exit set. It remains to decide in which of the K"-tiles the intersection of the tile with the exit set 
takes at least one half of the volume of the tile. 

To this end we use numerical simulations and the fact that for an autonomous system of ODEs 
x = f(x) with a solution x(t) the function x(—t) is a solution of autonomous system x = —f(x). For 
determining whether to include a fc-tile in the entry set of a successor state, we sample the tile and 
perform numerical simulations of the trajectories of system x = —f{x). If more than one half of the 
simulated trajectories go through the rectangle and the entry set of the original state, then the fc-tile is 
included in the entry set of successor state, otherwise the fc-tile is not included. 
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Weights of the transitions correspond to portions of the set of performed simulations that leave the 
rectangle to the respective neighbouring rectangles. Weight of the transition from the state to the so- 
called selfloop state with the same rectangle is determined as the portion of trajectories that do not leave 
the rectangle in given maximal time interval. 

Algorithm 2 Procedure getSuccessors 

Require: SS = (nj, ,T , J), K,MeN,He Rect(S), E 6 EntrySets(H) 
Ensure: Successors = Slices K ((H ,£)) 



1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 



if £ = then 

Successors <— {(H,&)} 
return Successors 
A <— set of M random points in E 
ExitPoints <- 
Stayslnside <— 
for all xq 6 A do 

simulate trajectory from xq until it leaves H through a point x\ or given time elapses 
if x\ exists then 

ExitPoints <— ExitPoints U{x[ } 
else 

Stayslnside <— Stayslnside +1 
for all FG Facets(H), ¥=HC\H' do 
if ExitPoints n F jt then 

Entry Tiles ^{Ze Tiles K (F) | Zn ExitPoints ^ 0} 
for all Z e Entry Tiles do 

B <— set of M random points in Z 
RealPointsCount <- 
for all yo 6 B do 

simulate trajectory from yo until it leaves H through a point yi or given time elapses 
if yi 6 E then 

RealPointsCount <— RealPointsCount + 1 
if RealPointsCount < f then 
Entry Tiles <— Entry Tiles \ {Z} 
ifEntryTiles^0then 

Successors <— Successors U(H' ,EntryTiles) 

Weighty, E)][(H,EntryTile S )] <- l ExitP °^ tsnF l 

return Successors 



Performing the backward simulations (lines 16-24 of Algorithm[2]) can be switched off. The resulting 



transition system differs from the QDAA in the entry sets, that can be larger. Difference of the outputs 
can be seen on Figure[5] The algorithm with backward simulations computes the QDAA and for (k — > °°) 
approaches the real behaviour of the solutions of dynamics ODE system. On the other hand the algorithm 
without backward simulations overapproximates the entry sets, therefore the transitions are included 
even if the entry set of a state is smaller than half of one K"-tile. Both options still lead to automatons 
with reachable states whose rectangles are included in the set of reachable rectangles of the rectangular 
abstraction with the same initial rectangles. 

The worst case complexity of the algorithms follows. There are at most k n rectangles in the phase 
space of the biochemical system, where k is the maximal number of thresholds on one variable. The 
maximal number of states of QDAA of the form (H,E) for a fixed rectangle H is 2n ■ (2 K " — l), where 
n is the dimension of the biochemical system. For the average numbers of visited different states of 
QDAA with the same rectangle encountered while analysing our evaluation models see the line labeled 
p in Table [T] Complexity of the computation of successors of a given state depends on the dimension 
of the system, the K parameter and on the number of simulations M used per one tile. In the worst case 
when all the tiles are examined (either as a part of entry set or potential exit set) there are In ■ K*" -1 -M 
simulations. 
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Visualization of the state space of QDAA involves highlighting the borders of the rectangles H such 
that there is at least one state (H,E) visited during the computation. The intensity of the fill colour of 
a rectangle H is calculated proportional to the sum of weights of all possible paths from initial set J^c 
to the first appearance of states with H as the rectangle. The weight of a finite path is obtained as the 
product of weights of the subsequent transitions in the path. The sum is always between zero and one. 

4.1 Entering and leaving conditions 

The computation of successor states proceeds in two steps. 

First, the probabilities of potential successors are computed for all fc-tiles of E and summed up to get 
the probabilities of successors for the whole set E. For each fc-tile of E several numerical simulations 
of a solution of system ([I]), with initial point x(0) placed randomly in the tile, are performed. If the 
computed trajectory satisfies entering and leaving conditions (Definition |4.2| ) and leaves the box H or the 
maximal time interval elapses, the number of trajectories leaving H through the particular facet (resp. 
the number of trajectories assumed to stay forever in H and leading to the transition {H,E} — > (//,0)) is 
increased. 

Second step of the algorithm takes into account the rectangles H' tx H with nonzero probability of 
transition {H,E} — > (H',yet unknown E') computed in the first step and determines the entry sets E' of 
this successors. Denote F =H' C\H. For every ff-tile of F the algorithm decides if the tile is a subset of 
E'. 



Replacing Entry(F,H) with more easily computed set Entry' (F,H) (defined in Definition 4.2 1 of 



points satisfying the entering condition in the definition of QDAA does not lead to a nonzero difference 



in the values of transitions probabilities due to Theorem 4. 1 



Definition 4.1 Let F G Facetst(H) be a facet of the form F = [ai,Z?i] x . . . {c} ... x [a n ,b n ]. Define the 
normal vector v# (F) to the facet F with respect to the rectangle H as the vector —eifor c = a, and eifor 
c = h[, where ei denotes the ith vector of the standard basis ofW (i.e. the vector orthogonal to F and 
pointing outside from H). 

Definition 4.2 Consider rectangle H,F,F' G Facets(H), a mutiaffine vector field f, a solution x(t) of 
the system HI and r > satisfying x(0) G F,x(r) G F',\/t G (0,r) : x(t) G H. We say that x(t) satisfies 
the entering (resp. leaving) condition with respect to H and f, if f(x(0)) • Vh(F) < 0) (resp. f(x(r)) ■ 
v H (F') > 0). 

Define Entry 1 \F,H) = {x G F\v H (F)-f(x) < 0}. 



Our next step is to show (in Theorem 4.1 ) that the focal set of Entry' (F,H) \ Entry(F,H) is of measure 
zero, thus the set difference is insignificant for our volume based notions. 

Now we will extend our definition of Focal(H,E,F') set of points from which the trajectories leave 



the box through a facet F' (Definition 3.6 1 to Focal(H,E,X) for an arbitrary subset X of facet F'. 



Definition 4.3 Let E G Entry Sets K (H),E / 0,F' G Facets(H). Let X C F'. Define Focal(H,E,X) the 
set of all points yo G E such that there exist £,c,e' > and a solution y(t) of the system nil satisfying 
y(t) ^H fort G (0,c),y(t) G Inter(H) for t G (0,£), y(c) G X, andy(t) G Inter (H') for t G (c,c + e'). 

Theorem 4.1 Let E G Entry Sets K (H),E / <d,F' G Facets(H). Let V+ = {x G ExitSet(H ,E ,F')\i] FI ■ 
f(x) > 0} and let Vq = ExitSet(H,E,F') \ V+. Then X*, (Focal(H,E,Vo)) = 0, where n' denotes nfor 
E = H, and (n — 1 ) otherwise. 
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Proof First, we have to observe that trajectories from such initial points have to leave the rectangle 
through a set of measure zero (in fact it is a zero set of a multi-affine polynomial). 



Then the proof is simillar to proof of Lemma 3.5 □ 



Remark 4.1 Theorem4.1 implies that replacing Entry (F,H) with Entry' (F,H) in the definition of QDAA 
(recall that these can be only entry sets of the successor states, not of the initial states (H,H)), does not 
lead to a nonzero difference in the values of transitions probabilites. 

Remark 4.2 Deciding if a point x is an element of Entry' (F, H) is straightforward (compared to checking 



the e-condition of Definition 3.2), and Entry' (F,H) C Entry(F,H). 



Remark 4.3 Moreover, there is the following symmetry property. Let F,F' G Facets(H),xo G F,xi G F'. 
There is a solution y(t) of system (T7I) satisfying y(0) = xo, 3?i : y(t\ ) = x\,f(y(0)) • v#(F) < 0,f(y(t\)) ■ 
Vh(F') > and y(t) G H for t G (0,?i), if and only if there is a solution x(t) of system x(t) = —f(x) 
satisfying x(0) = x u 3ti :x(t]) =x ,-f(x(0)) ■ V H (F') < 0,-/(x(?i)) • V#(F) > and x(t) G H for t G 
(CVi). 



5 Evaluation and Case Study 

In this section the state spaces of several biological models (of dimensions two, four and seven) are 
explored. Using our prototype implementation of the algorithms from Section [4] implemented in C++, 
we evaluate our approach on two exemplary biochemical systems. Additionally, we provide a case study 
held on a biochemical pathway studied in E. coli and compare the reachability results of the case study 
and one of the smaller models with results obtained using the rectangular abstraction approach. 

Before we proceed with the models, let us introduce several terms useful for the evaluation. For a 
biochemical system g$ = (n,f, £?, J?c) we denote M(^c) ^ Rect(^) the set of all rectangles reachable 
from initial set J?c- For each H G Rect(^) we denote mem(H) the subset of M(-J?c) consisting of 
all states reachable from the initial set with H as rectangle, the so-called memory of the rectangle H, 
mem(H) = {(R,E) G M(J?c) \ R = H}. Further we denote p the average number of memory states 
(cardinality of mem(H) averaged over all H G &(j?c))- The number of QDAA states representing the 
memory of a rectangle is in the worst case equal to the number of all its possible entry sets. However, 
the actual values of p in our examples are much smaller (see Table [T}. 

Let us focus on the effect of parameter K on cardinality of M(^c) and on p. Expected behaviour of 
the approximation is the following. Every facet is divided into K n ~ l tiles. A tile is included in the entry 
set E of some reachable state (H,E) if the focal subset Focal(H,E,A) fills at least half of the volume of 
the tile. For higher values of K, the set Tiles* (Focal (H,E, A)) better approximates the set Focal(H,E,A) 
because of the higher fc-grid resolution. Thus with increasing K, the quantitative information denoting 
the probability of reaching states in M(^c) can be computed more precisely. We demonstrate that on 
models examined below. 



5.1 Oscillatory model 

First, we consider a 2-dimensional model which is a variant of Lotka-Volterra model with oscillatory 
behaviour. The oscillatory model has the form of the following multi-affine system: 
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g=5-X-l-X.F 



%=0A-X-Y-5A-Y 



We consider the following partition and initial conditions for this model: 

Zi = {»|i€(0,30)CN } 
r y = {j|JG(0,12)CN } 

^c-X£ (20,21),7e(5,6) 

Results achieved on our implementation are presented in Tableland visualized in Figure [5] Black 
rectangles denote the initial set. 

5.2 Enzyme kinetics 

Similarly, we examined a 4-dimensional model of basic enzyme kinetics based on the following set of 
reactions: 

S + E ^>ES 
ES %S + E 
ES %P + E 
The corresponding multi-affine ODE model considered in the paper is the following: 

f = -0M-S-E + l-ES 
*§ = l-ES-O.Ol-S-E + 1-ES 
^§f = -\-ES + 0.0l-E-S-\-ES 
f = l-ES 

We consider the following partition and initial conditions for this model: 

T s = {0.01,5, 10, 15,25,50,60,85,95, 100} 
T E = {0.01,5, 10, 15,25,50,60,85,95, 100} 
T ES = {0.01,5,10,15,25,50,60,85,95,100} 
T P = {0.01,5, 10, 15,25,50,60,85,95, 100} 

J c : Se (25, 50), E £ (95, 100), ES e (0.01, 5), P 6(0.01,10) 

Projection of the approximated phase space to the enzyme/substrate plane is shown in Figure [6] For 
both the oscillatory model and the enzyme kinetics model full version of Algorithm [2] (with backward 
simulations) was used. 
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Figure 5: Reachability in oscillatory model and comparison with numerical simulation, first two figures 
were obtained using the full version of Algorithm |2| the third one with lines 16-24 omitted. For com- 
parison: using the rectangular abstraction transition system on this biochemical model, the whole phase 
space [0,30] x [0, 12] is reachable from the same inital conditions. 





Oscillatory 


Enzyme 


K 


4 


8 


16 


32 


64 


128 


4 


5 


6 


7 


\M{J?c)\ 
P 


52 
1.63 


46 

2.2 


40 

3.78 


39 
2.9 


37 
4.57 


35 
6 


76 
4.36 


104 
10.76 


123 
16.8 


166 

53.6 



Table 1 : Results for the two models and several different settings of the discretization parameter k. 




0.01 5 10 15 25 50 



0.01 5 10 15 25 50 



K = 4 



Numerical simulation 



Figure 6: Enzyme kinetics model - projection of the reachable set to the enzyme/substrate plane and 
comparison with numerical simulation. 



L.Brimetal. 21 



5.3 Case Study on E.Coli Ammonium Assimilation Model 

We consider a model specifying the ammonium transport from the external environment into cells of 
E. Coli [26]. The model describes the ammonium transport process that takes effect at very low ex- 
ternal ammonium concentrations. In such conditions, the transport process complements the deficient 
ammonium diffusion. The process is driven by a membrane-located ammonium transport protein AmtB 
that binds external ammonium cations NH 4 ex and uses their electrical potential to conduct NH 3 into the 
cytoplasm. In Figure [7] biochemical reactions of this model and the scheme of the transport channel are 
shown (left and middle). 

The level of pH and external ammonium concentration are considered constant. 

The system of differential equations: 

^^1 = -k\ [AmtB] [NH 4 ex] + fa [AmtB : NH 4 ] + fa [AmtB : NH 3 ] 

d[ Am tB:NH 3 ] = ^ y AmB . NH ^ _ ^ y AmtB . NH ^ 

d[AmtB:NH 4 ] _ ^ ^ AmtB ^ { NH4ex \ - fa [AmtB : NH 4 ] - fa [AmtB : NH 4 ] 
d[N " 3in] = fa [AmtB : NH 3 ] - fa [NH 3 in] [H in ] + fa [NH 4 in] + fa [NH 3 ex] 
d -^§^ = fa [NH 3 in] [H in ] - fa [NH 4 in] - fa [NH 4 in] 

Constant species: NH 3 ex,NH 4 ex,Hi„,H ex . 
Initial conditions and threshold numbers: 

T NH3ex = {0,28-10- 9 ,29-10- 9 ,l-10- 5 } 
T NHA ex - {0,49-10- 7 ,5-10- 6 ,l-10- 5 } 

10~ 12 , 1 • 10 10 ,5 • 10~ 6 ,9.9 • 10~ 6 , 1 • lCT 5 } 

10 7 ,1-10 5 } 

10 7 ,1-10 5 } 

10~ 8 , 1 • io~ 7 , 1 • 10 6 ,11 • 10~ 7 , 1 • 10~ 5 , 1 • io~ 4 , 1 • 10~ 3 } 

10~ 8 , 1 • 10~ 7 ,2 • 10 6 ,2.1 • 10~ 6 , 1 • 10~ 6 , 1 • 10~ 5 , 1 • 10~ 4 , 1 • 10~ 3 } 



TkmtB = {0, 1 

TAmtB.NHi = {0, 1 

TAmtB:NH 4 = {0, 1 

TNHiin = {0, 1 

TNH 4 in = {0, 1 



Jc- /V#3£*:e(28-10- 9 ,29-10- 9 ), 
NH 4 exe (49-10- 7 ,5-10- 6 ), 
AmtBe (0,M0~ 5 ), 
AmtB:NH 3 £ (0,1-10~ 5 ), 
AmtB:NH 4 e (0,1-10~ 5 ), 
NH 3 ine (1 • 10 6 ,11 • 10 7 ), 
NH 4 in£ (2-10- 6 ,21-10- 7 ). 

The upper bounds on concentrations of NH 3 in and NH 4 in considering the biological system with 
given initial conditions were estimated as 1.1 • 10~ 6 (NH 3 in does not exceed the initial concentration) 
and 5.4 • 10~ 4 by the rectangular abstraction (overapproximation). 

Reachable intervals using Algorithm ^without the backward simulations were [10~ 8 , 1.1 • 10~ 6 ] for 
NH 3 in (NH 3 in does not exceed the initial concentration), and [10 , 10 ] for NH 4 in. This results are in 
agreement with simulated data and in the case of the concentration of NH 4 in the QDAA results are by 
one order closer to numerical simulations than the rectangular abstraction results as can be seen in the 
right part of Figure [JJ 
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AmtB+NHnex £--4 AmtB : NH t ki=5- 10 s , k 2 = 5 ■ 10 3 

AmtB : NH A % AmtB :N Hi +H lx fe =50 

AmtB -.NH^h AmtB + NH 3 in fc, = 50 

OTJjfri 4 fcj = 80 



le-6 



NH d in+H ta ^-%NH 4 in k 6 = 1 ■ 10 15 ,jt 7 = 5.62- 10 5 /^\ le-7 I 

NH 3 in£-$NH 3 a h = h = 1.4- 10* /""' V \ le " 10 le-6 le-3 1 le+2 le+6 le+8 

«, / \^ time(s) 

Figure 7: Ammonium transport model (left). Simulations of the ammonium assimilation model from 20 
randomly sampled points in J 2 ^ projected on the concentration of NH^in, blue lines represent bounds on 
this concentration found by the QDAA - two subsequent thresholds 10~ 6 , 10~ 5 (right). 

6 Conclusion 

We have presented a new theoretical method for finite discrete approximation of autonomous contin- 
uous systems equipped with a measure that indirectly quantifies correspondence of the approximated 
behaviour with the original continuous behaviour. We have provided a computational technique which 
we implemented in a prototype software. We have examined the implementation on small dimensional 
models which showed satisfactory results for computing reachability. 

The method can be either used as a parameterized simulation technique or employed with rectangular 
abstraction to quantify the extent of spurious counterexamples. Thus the method can improve the current 
possibilities of analysis based on model checking techniques. We leave for future work integration of 
this method into the software for model checking of biochemical dynamical systems [13]. 

At the theoretical side, we leave for future work precise clarification of our method wrt the rectangu- 
lar abstraction. From the computational viewpoint, we aim to develop a parallel reachability algorithm 
that would make the method scalable and applicable to systems of larger dimensions. 
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